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ABSTRACT 


A model that produces three-dimensional ray trajectories in an optical fiber is 
derived through the use of a simple correspondence rule. The analysis and simulations 
presented will be in terms of dimensionless parameters. The curves prepared are 
presented for allowed radial and azimuthal mode numbers and are presented in terms of 
dimensions which are normalized by the core radius. The model presented will be shown 
to be in agreement with the standard Eikonal ray theory and will serve as a link between 


the ray trajectories and the mode numbers of the WKB waveguide solutions of an optical 


fiber. 
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I. INTRODUCTION 


A. REVIEW OF PRIOR WORK 

An extensive body of literature exists for both wave and 
ray characterization of optical fibers. It is generally 
accepted that exact waveguide solutions exist for a select 
number of specific index profiles [Ref. 1, 2]. The exact 
waveguide solutions are frequently derived for step index 
fibers in intermediate level optical fiber texts. For example, 
Cherin’s Optical Fiber Communications [Ref. 3] provides a 
fairly clear waveguide analysis of the step index fiber. In 
the general case of graded index fibers the problem is not 
analytically tractable without approximation [Ref. 1]. The 
most common technique for the analysis of graded index fibers 
is based on the WKB approximation (Ref. 4, 5, 6, 7]. It is 
b.ssible to show that the first-order WKB approximation is 
equivalent to the results produced from a simple ray model 
(Ref. 8]. 

Despite the mathematical importance of having an exact 
solution for the waveguide fibers in an optical fiber, the 
propagation characteristics can be difficult to visualize. It 
is not surprising that an alternative approach to treating 


optical fibers, based on the Eikonal ray approach [Ref. 9, 10, 


11, 12), which in turn is based on Fermat’s extremum principle 





(Ref. 10, 11, 12, 13, 14], has proven to be extremely useful 
in explaining effects which take place in optical fibers. For 
example, application of such techniques to the prediction of 
bending losses in optical fibers has appeared in the 
literature (Ref. 13, 15]. In addition, achromatic modal 
dispersion effects have been predicted with ray models [Ref. 
16). It has been shown that the exact step index radical 
function has a first inflection point at the inner ray caustic 
or turning point [Ref. 8]. To a significant degree the 
association between the WKB wave solutions, which depends on 
azimuthal and radial mode numbers, and ray models has not been 
fully explored. For example, the ray trajectories were 
typically derived for arbitrary initial conditions without 
regard for ray conditions dictated by invariant conditions 
associated with mode specific quantum numbers. 

Recently, an approach based on a simple correspondence 
between waves and rays has been proposed [Ref. 17] as a 
substitute for the first-order WKB approach of counting modes 
[Ref. 9]. This correspondence principle has been used to 
qualitatively demonstrate that there should be a direct 
association between the mode numbers in an optical fiber and 
the ray trajectories. The main focus of this thesis is to 


treat this association quantitatively. 


B. OBJECTIVES OF THE THESIS 

The correspondence principle [Ref. 17] leads to rules 
which predict both initial conditions and propagator 
conditions on ray trajectories associated with specific modes. 
This thesis will cover all background analysis leading to 


application of these rules to meridional rays and skew rays in 


a graded index fiber. Whenever necessary for visualization, 


full three-dimensional ray trajectories are provided. The 


analysis and program simulations are derived and implemented, 





respectively, in a normalized form with respect to the vacuum 
wavelength and core radius. As a result, the computer results 
are generalized to a broader class of problems. Specifically, 
the optical fiber rays are expressed in terms of the commonly 
used V-parameter, numerical aperture and core index. Appendix 


B has been dedicated to demonstrating that the correspondence 





principle, mentioned above, leads to propagator conditions 
which are precisely equivalent to rules obtained from the 
standard Eikonal model. Appendix A has outlined in flowchart 
fashion the conceptional development of the programs. Appendix 
E of the thesis manuscript has been scheduled for an IEEE 


publication [Ref. 19] in March 1993. 


C. OVERVIEW OF THE THESIS 
Chapter I is the introduction and covers both a review of 
relevant literature and a thesis overview. In the interest of 


keeping this thesis fairly self-contained, Chapter II provides 


a brief technical description of essential background 
preliminaries. 

In Chapter III, the basis for the specific problems 
examined in later chapters is provided. A conceptual flowchart 
shown in Appendix A summarizes the basic analysis which leads 
to mode specific trajectories. 

In Chapter IV, attention is given to meridional rays for 
the arbitrary oa profile. The specific meridional ray cases of 
step index fibers (a=) and parabolicity (a=2) are both 
analyzed in detail. Appendix D provides a useful mathematical 
analysis for the essential function F(a), a normalization 
integral used in this chapter. 

In Chapter V, the solutions for the step index skew rays 
are obtained. This chapter is a link between the general 
formulation of Chapter III and the more complete step index 
analysis contained in Appendix E. 

In Chapter VI, numerical solutions for skew rays with 
arbitrary @ profile is provided. Some useful variations of 
fiber optic parameters due to initial conditions are also 
demonstrated and explained here. 

In Chapter VII, quasi analytic solutions for the 
interesting case of skew rays with a=2 are presented. The well 
known fact of the orbital periodicity of this case is also 
proven here. 

Appendix B shows the exact agreement between the ray 


analysis presented and the Eikonal analysis. 


Appendix C provides a symbol table which contains notation 


for the most common symbols used throughout in the thesis. 





TI. BACKGROUND 


Starting from the Hemholtz equation, which is the phasor 
form of the linear wave equation for isotropic, source-free, 


homogeneous media, we have 

[V?+k? (x) ] B=0, (2-1) 
where k(r) is the wave number and E is the electric field. It 
is noted that this equation is only approximately valid ina 


graded index fiber. Consistent with the form given in (2-1), 


we assume that 
k (r) =n(r) k=k,6,+k,6,+k,é,, (2-2) 


where k is the vacuum wave number, n(¥) is the index of 
refraction, (r, ¢, z) are the cylindrical coordinates needed 


to describe points in the fiber and (é,, 6,, 6,) are the unit 


vectors of the coordinate system. 
We list here the basic assumptions that n(7r) satisfies 


1. The index depends on the distance from the symmetry axis 
of the fiber (viz. n(r)=n(r)). 


2. The index is a decreasing function for rsa. 
3. The index is constant for rea. 
These restrictions can be represented in Figure 1 where 


K(r)=kn(r). 














Figure 1: Restrictions for k(r). 


Following a standard procedure [Ref. 1] in waveguide 
analysis, we assume that the general phasor for the z 
component of the electric field in the cylindrical optical 


fiber can be represented as 


E,(r) =F,(r) ee, (2-3) 


where F,(r) is the radial solution to (2-1), 8 is the phase 
constant of the waveguide, and v is the azimuthal mode number, 
which is forced to be an integer due to periodic boundary 


conditions [Ref. 1]. 





In the case of a step index fiber, an application of 
Maxwell's equations leads to a transcendental equation for the 
phase constant. The equation has a finite number of discrete 
solutions for each integer »y [Ref. 1]. The obtained convention 


for ordering these solutions is 


B,>8,..->B,>B,_. (2-4) 
where the second subscript m is known as the radial mode 
number. This convention is also applied to the approximate WKB 
(Ref. 7] solutions in inhomogeneous optical fibers. 


A simple correspondence rule for plane waves [Ref. 17] 
v=-jk, (2-5) 


can be applied to the ¢ and z dependencies of (2-3). The grad 


operator is given by [Ref. 18] 
726 foe 2 (2-6) 
It now follows from (2-6) and (2-3) that 
(VE) =J2E,, (2-7) 
which from (2-5) leads to 
: (2-8) 


Similarly, 
j (VE,) ,=-JBE,, (2-9) 


and again from (2-5) 





k,=B. (2-10) 


It is noted from (2-3) that the correspondence rule applies 
only to the ¢ and z components. 

Using the relation between the propagation vector and its 
components, it follows from (2-2), (2-8) and (2-10) that the 


radial component of the propagation vector can be expressed as 
ree 222 -B?-(2)?. (2-11) 
, r 


Both + and - roots are required for the development of a full 
cycle of a trajectory. Nevertheless, symmetry conditions in 
the trajectory will allow the analysis to be based on the + 
root only. Furthermore, for k, to correspond to a bounded 
orbit, it must be real [Ref. 1, 8). Without loss in 
generality, we can continue with only the + form for k,. 
According to (2-4) and (2-11), the radial component of the 
propagation vector will increase with increasing m (see Figure 
2). The crosshatched domain of Figure 3 represents the 
condition that k, is real. Therefore, terms under the radical 


of (2-11) satisfy 
k?=k? (r) -6?-(2)?20. 
r (r) B (—) 2 (2-12) 
A condition of phase synchronization [Ref. 1] on the 


radial part of the wave vector results in the approximate 


expression 





\ 


. KO OOAK s = 
increasing | Bel \ increasing m 


kK: 


decreasing 8 








Figure 2: Radial component of propagation vector k, versus m. 
m2f*| k? (r) -6?- (2)? ar, (2-13) 
T Jr, r 


where the radial mode number m must take on positive integer 
values. The integral limits r, and r,, are the two values for 
which the integrand of (2-13) is zero. These are generally 
referred to as the turning points as represented in Figure 3. 


We note that r, and r, satisfy 
kp =k? (x) -8?-(—)?=0. (2-14) 


More precisely, the relation (2-13) becomes exact after the 


1 


substitutions mom = and pops (Ref. 5]. In the WKB analysis 


of optical fibers, these offsets are usually ignored. 














Figure 3: Condition for bounded trajectory. 


A significant point relates Figure 3 to the radial mode 
number given by (2-13). The value of m is qualitatively 
proportional to the cross-hatched area. The general procedure 
to find the maximum m is to allow f? to be equal to k,’. If Bf? 
falls below k,’, then the additional modes created are no 
longer trapped. Many authors have explored the interpretation 
that the optical fiber trajectory is equivalent to the 
trajectory of a particle in a potential well [Ref. 8]. 

Figure 3 correctly suggests that there is a maximum p for 


which the cross-hatched region vanishes. Clearly there is an 


11 


upper limit on the value of v which occurs for the condition 
r,=r,. The general relation between the turning points and »p 


can be obtained from (2-14) as 


y=xVk?(r) -6? , iar) 


which is represented in Figure 4. It is apparent from Figure 
3 that the condition #=kn, will allow the maximum range in p 
and the corresponding maximum in the peak on Figure 4. In 
calculating »,,,, v(r,) will in general be evaluated under the 
stated limiting condition fS=kn,. At v=0, the solutions for 


turning points are easily read off in Figure 4 as 0 and 


BUE Ys Since 8 <n,cos6, the second condition leads to the 


condition for total internal reflection 

n(r,) =n,cosé. (2-16) 
For the step index fiber, r,=a, which leads to the well known 
relation 


cos@=_2, (2-17) 
1 


Equations (2-8), (2-10), (2-11), (2-12) and (2-15) become 
the basis for the three-dimensional ray modeling that will be 
the focus of the discussion that follows. It is noted that the 
components (k,, ky, k,) will describe the direction of the ray 


and that coordinates (r, ¢, z) locate the ray. For purposes of 


12 
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Imr. GENERAL CHARACTERIZATION OF THE THREE-DIMENSIONAL RAY 


MODEL 


A. INTRODUCTORY REMARKS 

This chapter will provide a basis for the specific 
categories of problems examined in later chapters (e.g., 
meridional, skew step index, and skew graded index). A general 
approach to the goal of obtaining mode-connected ray 


trajectories will be presented at the end of the chapter. 


B. GEOMETRY OF THE PROBLEM 
With reference to Figure 5, the orientation of the wave 


vector is defined in terms of @ and ~. It follows that 


k,=k(r)sin@ sing, (3-1) 
k,=k(r)cos6, (3-2) 

and 
k,=k(r)sin@ cosé. (3-3) 


The incremental line variation in the ray trajectory can be 
defined as 

dl=dz é@,.+dr 6 +r dg 6,. (3-4) 
The propagation vector, defined by the cylindrical components 


(3-1), (3-2) and (3-3), is perpendicular to the wavefronts of 


constant phase and therefore tangent to the ray trajectories. 


14 





Figure 5: Geometry of the problem. 


15 








Therefore, the incremental components dd, dr, and dz are 


tangential to (2-8), (2-10) and (2-11), respectively, and it 
foliows 
dz=dl cos@, (3-5) 
dr=dl sin@é cosé (3-6) 
and 
r dod=dl1 siné siné. (S27) 


C. MODELING ANALYSIS 

1. Computation of the turning points 

Without any loss in generality, we assume that the index 
profile can be expressed in terms of a standard power profile, 


which is matched to the edge conditions, 

n(0) =n, (328) 
and 

n(a) =n, (3-9) 
and a perturbation f(r) as shown below 


k'n?(r) =k? (n} (x) +NA7E(4r) 1, (3-10) 


where 


£(0) =f(a) =0. (3-11) 


The numerical aperture of the fiber is given by [Ref. 1] 
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NA= n-n} (3-12) 
and the power profile satisfies 
ny=mj [1-2A(=)], (3-13) 


in which the index difference for this formula is defined as 





2 2 
hee (3-14) 
2n; 
It can be easily checked that at r=a, n,=n,. 
Substitution of (3-12) into (2-12) leads to 
w LE meen (bre E ioe) 
T dr, V r 
where the V-parameter is given by [Ref. 1) 
V=kNAa, (3-16) 
the normalized phase constant, ¢€(8), is defined by 
2 B 2 
n,- (=) . 
€(B) = k (3-17) 
NA? 
and 
ro (3-18) 


a 


defines the radial normalization. 


17 





In order to further analyze the problem, the perturbation 
of the index power will be dropped. After taking f(r)=0, the 


condition on the turning points will satisfy 
re?—¢ (B) r?4(—)?=0, (3-19) 


which is a normalized version of condition (2-14). Figure 6 is 
plotted for a@=2, V=100, v=10 and ¢(f8) as stated on the Figure. 
It is noted that as ¢(f8) increases, the zero crossing points 
(or turning points) should move closer together. Eventually 
the turning points will coalesce. 

2. Computation of the maximum radial mode number 

As described in Chapter II, the maximum radial mode number 
is obtained by setting B=k,, or according to (3-17), ¢€(8)=1. 


Thus, (3-15) leads to 
ai V (1 S27 By2,5e2 dr (3-20) 
Max vel (+) +x ae 


where the turning points which correspond to the maximum m 


will then satisfy 
Fr 24 (= )2=0, (3-21) 


3. Computation of the maximum azimuthal mode number 
The first step in obtaining the maximum azimuthal mode 
number »v,,, is to find v(r,) as represented in Figure 4. There 


is a peak in the function pv(r) already defined by (2-15). The 


second step is to evaluate the expression for the limiting 





nornalized radius 





Figure 6: Left hand side of Equation (3-19). 


condition #=kn,. 


Combining the power profile (3-13) with (2-15) leads to 


yt=rAk2n? [1-20 ( = -p, (3-22) 
By evaluating 
d 2) = 3- 
7 oe, 0, (3-23) 


it can be shown that 


(arbitrary v), ree) 





which for the limit condition B=kn,, leads to 





1 
Ste ee 2 (3-25) 
Zz (sa) (v=P,,.)» 
where z. is the normalized radius extremum given by 
rele (3-26) 
pa 


Finally, from substitution of (3-25) into (3-24), we have 








1 
= 2 )\% a (3-27) 
Paar V( oT) jane 


which shows that »,,, depends on the qa profile and the V- 
parameter. 
4. Computation of the evolution equations 


Through combinations of (3-2), (3-3), (3-5) and (3-6) we 


have that 
dz_dz_*,_ 6 _ctné (3-28) 
dr dr k, k, cost’ 
where 
2-2 (3-29) 
a 
is the normalized z component of the ray. 
It is easy to demonstrate from (2-11), (2-12) and (3-15) 


that 
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eo ¥| care etn) 2), (3-30) 
r v 


and from (3-16) and (3-17) that 


B yn; -NA?e (B) (3-32) 


zi a NA : 


Substitution of (3-30) and (3-31) into (3-28) leads to 


n?-NA7e (8) E 


NA — —_— 
[« (B) r?-(7)?- (x) 2+ 


ala 


f (3-32) 


which is the evolution equation for the normalized z component 
of the ray. 
Working similarly for the normalized ¢ component of the 


ray, combination of (3-1), (3-3), (3-6) and (3-7) gives 


ee (4233) 
| ean -BP- (=)? 


Kis 


and after substitution of (3-30) leads to 


va 1 

VE arena (3-34) 
e(B)r ~ Cy) -r 

which is the evolution equation for the normalized ¢ component 

of the ray. Equations (3-32) and (3-34), which describe the 

Yay trajectories, are shown to be in agreement with the 


Eikonal solution [Ref. 9] in Appendix B. Unlike the Eikonal 


21 





evolution equations, (3-32) and (3-34) are expressed in terms 
of the azimuthal mode number and the phase constant associated 
with the radial mode number (2-13). 


The procedure describing the general approach applied in 


subsequent chapters is summarized in Appendix A. 





IV. MERIDIONAL RAYS 


In this chapter attention is given to meridional rays. It 
is well known that modal dispersion can be approximated fairly 
well through consideration of the meridional rays only [Ref. 
1]. In the first section, the treatment for an arbitrary a 
parameter is presented. In the last section, some conclusions 


for the specific case of step index fibers are then presented. 


A. GRADED INDEX MERIDIONAL RAYS WITH ARBITRARY a 
Meridional rays are rays which pass through the z axis, 


and according to (2-8), have no rotational component, that is, 
v=00k,=0. (4-1) 


For a power profile given by (3-13), it follows from (4-1) 


that (3-15) and (3-20) simplify for this case to 
mv (Ve (B)re-re? GF (4-2) 
Tt Jo x 
and 
Man ( Ege (4-3) 
r 


respectively. It is noted that the lower limit of integration 
in (4-2) and (4-3) is zero since the associated integrands 


have a root at r=0. This agrees with the physical 


interpretation consistent with Figure 7, which represents the 














Figure 7: Physical interpretation for meridional rays. 


special case v=0 of Figure 3. Equation (4-2) can be rewritten 


as 


, 


a ae (4-4) 
: € (B) al 


meV VE AB) ial 
us 9 


and therefore, the second turning point will, from (4-4), 


satisfy 


1 
F,=(€(B)) *. (4-5) 


After the variable substitution 
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oS et (4-6) 


1 
> (€(B))@ 
we find 
Qea 
Sa - 
ioe VETO) Pia), (4-7) 
where 
F(a)={ 1-u* du. (4-8) 


Although a general solution to (4-8) for all @ does not appear 
possible, the integral F(a) is easily computed via numerical 
simulation (see Figure 8). The functions F,(a) and F,(q) 
defined in Table 1 are seen in Figure 8 to agree with F(a) for 
small and large a, respectively. The analytic properties of 
F(a) are discussed in detail in Appendix D. Table 1 


summarizing this discussion is provided here. 





ee ees eS eae 
AS AOS SOR” ARES 


ae ae 
ae See 0_(non- susie 








F(a) = (2a) / (202) 
alam /2_ | 















Table 1: Values for F(a). 


25 


3 
Q 
a 
a 
M4 
Vy 
> 
3 
hy 


Figure 8 





In order to evaluate m,,, from (4-3), €(8) is again set to 


one (see Chapter III Section B), in (4-7), leading to 


Myx = F(t) « (4-9) 


A representative meridional graded index ray trajectory is 
shown in Figure 9. 
Before deriving evolution equations for ray trajectories, 


we observe from (2-10) and (3-2) that 


B=k n(r) cosé@,, (r) =k n,cos6,,, 





Figure 9: Representative meridional ray trajectory. 


27 





where the subscript "1" is to indicate evaluation at n(r)=n, 
(see (3-8)), and subscript m is to indicate a discrete angle 
associated with a specific radial mode. A geometrical 
interpretation which relates the angle 6,, to m can be made. 
After substituting (4-10) into (3-17) we get 


28 
n,sin’é,,, 
NA? 


€ (B) = (4-11) 


’ 


which when substituted into (4-7) provides a relationship 


between 6,, and m as 


a 


_NA,amm 1 ) = 


i ao (oe (4-12) 
sind,,, a '-V Fla) 





, 


which is in agreement with Figure 2. 
For calculation of z(r), we see that (3-32) simplifies 


under condition (4-1) to 


dz_\n;+NA’e (B) r (4-13) 


ar NA le (B) r7- (r) ° 


After combinations of (3-17), (4-10), and (4-11), z(r) takes 


on the integral form 


ae _n,cos8,,, 77 dr 
n,sin6,, Jo Pe fo a (4-14) 
1-(2) 
tT, 


After use of (4-6), it is easily seen that 
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2(F) =F,ctn8,,(, Les, (4-15) 


¥1-u* 


where re[0,r,]. The parameter r, is defined in terms of m 


through (4-5), (4-11), and (4-12) which leads to 


H(t. (4-16) 
Again, as in the case of F(a), exact solutions for (4-14) and 
(4-15) can be found for a=1, a=2, and a=o. 

Because of condition (4-1), the azimuthal variation is 
zero. Therefore, there is no need to consider the ¢ variation 
described by (3-34). Appendix A summarizes the process of 
trajectory analysis for meridional rays also. 

Figure 10 illustrates meridional rays plotted using (4-14) 
under conditions a=6, NA=0.5, V=100 and n,=0.5. The curves are 
distinguished only by the radial mode numbers m=20 and m=2. As 
expected from Figure 2, m increases with k,, and from (2-4), 
8 decreases with m. Examination of Figure 6 confirms that the 
peak values of the trajectories increase with m. The values 
predicted from (4-16), r,=0.9114 for m=20 and r,=0.5125 for 


=2, agree with the curves shown. 


B. GRADED INDEX MERIDIONAL RAYS WITH ax2 
Equations (4-9), (4-12), (4-15) and (4-16) can be 
evaluated for a=2. From Table 1, (4-9), (4-12), and (4-16) 


become 
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Figure 10: Meridional ray trajectories for a=6. 








Mux =Int 2) , (4-17) 
j _ NA 4m (4-18) 
siné,,, ale 
and 
> .| 4m (4-19) 
af. 


respectively. And lastly, through use of the integral [Ref. 


18] 


1 
f au _ -arcsinu, (4-20) 
0 2 

l-u 


the trajectory equation (4-15) can be expressed as 


rer,sin(——* __). 


——— (4-21) 
r,ctné,,, 


Therefore, for a=2, the trajectory shape is exactly sinusoidal 


with a normalized z period 


A,,=27r,ctn (4-22) 


Im? 
defined in terms of the specific initial condition of the 
angle of the ray as it crosses the normalized z axis. From (4- 


18) and (4-19), the dependence of the normalized period on 


V,NA and m is given by 











yt): (4-25) 
NA 


Figure 11 is a plot of A, for typical values of X. 

Figure 12 shows two meridional rays plotted using (4-22) 
under conditions a=2, NA=0.5, V=100 and n,=1.5. The curves are 
only distinguished by the radial mode numbers m=20 and m=2. As 
expected from (4-19), the peak value of the trajectories are 
T,=0.894 for m=20 and r,=0.282 for m=2. Also, it can be seen 
from the X=9 curve of Figure 11, that the period for m=2 is 


approximately 19 and the period for m=20 is 18. 


C. STEP INDEX MERIDIONAL RAYS 
By definition, for step index meridional rays, in the 
limit as 
Ooo, (4-26) 
which when substituted into (3-27) gives 
Van =Int(V)} . (4-27) 


The peak turning point can be found by substitutiig (4-1) 


and (4-26) into (3-25) yielding 








£ 
+ 
© 
a 
£ 
5 
c 
w 
0 
oO 
£ 
Ke) 
© 
ee 
5 
£ 
N 
rs) 


(@) to + N 
N = eo = 


poiiad jopiosnuts 





Figure 11: A, versus m for meridional rays with a=2. 
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Figure 12: Meridional ray trajectories for a=2. 








Y,=7,=1. (4-28) 
In order to find the maximum radial mode number (4-1), (4- 
26) and (4-28) are substituted into (3-20) which gives 
V 
Myx=Int[—) . (4-29) 
The evolution equation for the z component can be derived 
by substituting (4-26) into (3-28) yielding 


Gz _otné. (4-30) 
dé 


The still unknown angle @ of (4-30) can be found from (2- 


12), which for this case takes the form 
= 5 [jamsinedr. (4-31) 


By integrating (4-31) and eliminating a and k through 


combinations of (3-12) and (3-16), we have 





-eiy-1, MNAT £ 
6=sin" ( a (4-32) 
Finally, by integration of (4-30), we get 
z(r)=retné, (4-33) 


where according to (4-28), re[0,1]. Equations (4-32) and (4- 
33) are the basis for the simulations described in the next 
paragraph. 

Figure 13 illustrates step index meridional rays plotted 


using (4-33) under conditions NA=0.5, V=100 and n,=1.5. The 
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Figure 13: Step index meridional ray trajectories. 
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curves are only distinguished by the radial mode values m=1 
and m=20. Although it is not apparent from Figure 13, as 


expected from (4-28), in both cases, r,=1. 


D. FINAL REMARKS 

In this chapter, the general equations of Chapter III are 
applied under conditions valid for meridional rays. An 
analytic form for the ray trajectories was established. The 
existence of a direct connection between the radial mode 
number m and the orientation angle @ has been demonstrated. 
This is a first sign that a link between Eikonal ray 
trajectories and mode numbers is possible. It is noted that 
although a three-dimensional perspective of meridional ray 
trajectories could have been demonstrated, the fact that the 
rotational component is zero (see (4-1)) made it unnecessary. 
In the chapters that follow, the analysis is applied to skew 


rays. 
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Vv. SKEW RAYS IN A STEP INDEX OPTICAL FIBER 


A. INTRODUCTORY REMARKS 

This Chapter presents a brief discussion for the analysis 
of skew rays in a step index fiber. It follows from (3-13) 
that the general results from Chapter III can be applied after 
taking the limit a. A more detailed description of this 
subject can be found in Appendix E. This appendix is a 
reproduction of a published work [Ref. 19] based on results 
obtained in the early stage of the thesis. One purpose of this 
chapter is to provide a clear link between the general 
formulation of Chapter III and the analytical approach for 
step index optical fibers described in Appendix E. Some 
additional results will also be presented. In this chapter the 
notation on the turning points will be modified to be in 
agreement with Appendix E. Specifically, the more suggestive 


labels of r,, and r,,, will be used for r, and r,, respectively. 


B. OVERALL ANALYSIS 
In this section the solutions in the step index case (in 
the limit as ao) are obtained. It follows from Figure 14 that 


Yma=1-. This can be confirmed from (3-25) 
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Figure 14: Physical interpretation for step index skew rays. 


I 
> lim, 2 ,4 (5-1) 
Faux“orom | 55a) “*2- 


It then follows from (3-15), after taking f(r)=0, that 
oe At: i Piso ae (5-2) 
m2 lear (5) dr, 

where the lower limit can be obtained from (3-19) 


Fo emieenid 


nin = —————— - (5-3) 
We (B) 


After substitution of (5-3) into (5-2), we find that 
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which can be directly integrated using [Ref. 18] 


J 2_ 2 
{Ay auevut-a -asec™! (=) : (5-5) 
which leads to 


mel (ghd meet, (5-6) 


ro. Do 


which apparently is not applicable for meridional rays (v=0) 
since r,, given by (5-3) for meridional rays is zero. This 
case has been covered in Chapter IV, Section C. 

As in Chapter III, Section B, the computation of m,, 
requires setting ¢€(f8)=1.0. It follows from (5-2) that for 


M=Magy 


i yp _ 


and the corresponding value of m,,, is 
My (9) = = (WV ¥? -ysec" (©) I, (5-8) 


obtained after substituting (5-7) into (5-6). 


From (3-27) the maximum azimuthal mode number is obtained 


Vnag= Vi (5-9) 


in the limit as qo. 
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This section is provided as a link between the general 
formulation of Chapter III and the more detailed step index 
analysis contained in Appendix E. It is noted that (5-3), (5- 
7) and (5-8) are easily shown to be in agreement with (E-20), 


(E-21) and (E-23), respectively. 


C. ADDITIONAL RESULTS 

The results and discussion appearing in Appendix E were 
presented at an IEEE conference [Ref. 19]. Due to the 
limitation of space, only a set of representative trajectory 
curves were presented at the conference. The parameters for 
other trajectory examples are summarized in Table 7 in 
Appendix E. In order to illustrate the application of the 
program under limiting conditions, two additional examples are 
presented. 

1. Example A. v=1<p,,,=150 

As seen in Figures 15 and 16, the rays are nearly 


meridional. The value for 1£,,,=0.0104, which for m=30 and y=1, 


can be checked by use of (5-6). It is noted that the ¢ 


dependence on r is given by 


g=sec!(_=). (5-10) 


Fnin 


The trajectory lines are absolutely straight as expected ina 
step index fiber. In Figure 15, a full three-dimensional view 


of the trajectory is provided. It is noted that the 
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Figure 15: Step index skew rays with p<p,,,. 
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Figure 16: Projection of Figure 15 onto a polar plane. 
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trajectories continue out to r=1, which can be seen from 
Figure 16. 

2. Example B. v=191=sy,,.=200 

Here »y was chosen to yield a value r,, that was much 
closer to one, i.e., a high rotational component. Again for 
M=1 and v=191, Tyin=0.9358 can be checked with (5-6) and Fiyure 


17. The corresponding three-dimensional plot is Figure 18. 


D. FINAL REMARKS 

In this chapter ray trajectories were established from the 
general analysis of Chapter III with the only assumption that 
aso. From Appendix E, we can see that the model relating the 
mode numbers m and pv with the orientation angles @ and é is 
very simple. Also, the three-dimensional perspective of ray 
trajectories shows that the model provides a good physical 


insight for a deeper understanding of optical fibers. 
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Figure 17: Polar projection for step index skew rays with 


VV max 
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Figure 18: Three-dimensional representation of Figure 17. 





VI. NUMERICAL SOLUTION FOR SKEW RAYS WITH ARBITRARY a 


A. INTRODUCTORY REMARKS 

The treatment for skew rays with arbitrary oa requires the 
full generality of the analysis presented in Chapter III. 
Meridional rays, which are covered in Chapter IV, will not be 


treated here. 


B. ANALYSIS 

For the power profile given by (3-13), maximum azimuthal 
and radial mode numbers can be predicted from (3-27) and (3- 
20), respectively. Integral limits for (3-20) will satisfy (3- 
21). 

In order to characterize the trajectories given by (3-32) 
and (3-34) associated with the »v and m modes, the normalized 
phase constant given by (3-17) and the turning points need to 
be accurately determined. In the most general case the problem 
does not have an analytic solution. In fact, the condition on 


the turning points r, and r,, the normalized phase constant 


€(8), and the radial mode number given by (3-15) need to be 

satisfied in a self consistent manner. The conceptual 

flowchart for this process is provided in Figure 19. 
Finally, evolution equations can be predicted from the 


integral form of (3-32) and (3-34), which become 
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Given msm, 
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€, €[0,1] 





Find turning points (3-19) 
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Compare m with m 
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Figure 19: Calculation of e¢(8) and turning points for skew 
rays with arbitrary a. 
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/n?-Na? = an 
rae n, -NA’e (8) f rdr 


NA if ooo (6-1) 
[eur (5)? - (rx)? 
and 
dr 
$(r) = =F (6-2) 
€ (B) r?-(— a 


respectively, where re[r,,r,]. 


C. EXAMPLES 

1. Ray trajectories 

Figure 20 shows a three-dimensional plot of a numerically 
calculated skew ray trajectory from (6-1) and (6-2) using as6. 


The dense vertical lines are drawn at r-=r,, and the low 
density lines are drawn at r=r,. From Figure 20, Az is defined 


as the axial distance between ray turning points r, and r. 
Because of the three-dimensional perspective, the normalized 
increment in Az only appears to change size. Figure 21 shows 


the corresponding axial view. From Figure 21, A¢ is defined as 


A¢= saeocant FAS, (6-3) 
Yr, 


Figure 22 shows a ray trajectory with conditions given on 


the Figure with p~p,,,. As can be seen from Figure 23, which is 








Figure 20: Ray trajectory for conditions as stated. 
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the axial view of Figure 22, for this case, r, is close to Yr, 
as expected from Figure 3 for v close to »,,,. 

For the ray trajectories shown in Figure 24 and Figure 25, 
a very high value of a@ was chosen. As can be seen, the 
trajectories are indistinguishable from straight lines, as 
expected for rays in a step index fiber (as discussed in 
Chapter V for step index case in the limit as qo). 

The equations (6-1) and (6-2) were also numerically 
integrated for the case a=2. The three-dimensional perspective 
and axial views are shown on Figures 26 and 27, respectively. 
It is noted that in Figure 27, it is clear that the azimuthal 
orbit is closed. These results will be reexamined and 
explained in the next chapter. 

2. Variation of fiber parameters 

Table 2, Table 3 and Table 4 were derived from repeated 
applications of the optical conditions V=100, NA=0.5, n,=1.5 
and qa as stated on the tables. The corresponding three- 
dimensional representations are seen in Figures 26, 20 and 22, 
respectively. The trends relating A¢ and Az are readily 
explained. Because all of these tables have the same kind of 
variation of parameters, discussion will be concentrated on 
Table 2. In the lower half of the table, mis held constant at 
the value 10 while v varies between 4 and 16. As v increases, 
the angle Ad increases with the rotational component, k,, as 


expected. It is clear from (2-13) (or in normalized form from 
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Figure 24: Graded index ray trajectory with very high a. 
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Figure 25: Projection of Figure 24 onto a polar plane. 
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Figure 26: Ray trajectory for a=2. 
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Figure 27: Axial view of Figure 26. 
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om fe ee me | || treat 
}2 | 20 | so | 2 | 0.327 | 0.620 | 9.170 | 0.904 | 
js | 20 {| so {| 4 | 0.209 | 0.609 | 9.126 | 0.794 | 
je | 20 | so | 6 | 0.265 | 0.754 | 9.003 | 0.675 _| 
j20| 20 | so | 6 | 0.26 | 0.822 | 9.039 | 0.580 | 
p20 | 4 | so | 23 | 0.057 | 0.690 | 9.270 | 0.264 | 
jio| 6 | so | 22 | 0.208 | 0.740 

19 5 . 

17 





Jao | 12 | so | a9 | c-as2 | 0.785 | 9.083 | 0.382 | 


Table 2: Variations of fiber parameters for a=2. 

(3-15) and (3-17) that with m fixed any increase in vy must be 
compensated by a decrease in 8. The interpretation of as the 
z component of the propagation vector K(r) (see Figure 2) then 


clarifies the slight decrease in Az. In the upper half of the 


table, v is held fixed while m varies between 2 and 10. 





Table 3: Variations of fiber parameters for a=10000. 
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fe Pg tr ee | | ae) ap re | 
| 6s | 5 | 26 | 5 | 0.239 | 0.705 | 24.075 | 0.396 
| 68 | 30 | a5 | 5 | 0.425 | 0.850 | @.229 | 1.030 _| 


| 6s | 45 | 2 | 5 | o.si2 | 0.902 | 6.745 | 2.276 | 
| 6s | 5 | 26 | 25 | 0.052 | 0.906 | 7.083 | 0.105 
| 6a | 5 | 26 | a5 | 0.073 | 0.279 | 9.128 | 0.166 
se | os | 26 | 5 | 0-239 | 0.705 [14.075 | 0.396 





Table 4: Variations of fiber parameters for a=6. 


According to Figure 2, it is clear that as m increases, 8 
decreases, and as such, Az should decrease. This pattern is 
again evident from the table. Referring back to the relation 
between B and ¢(8), given by (3-17), and the discussion of the 
dependence of turning points on ¢(8) in Chapter III, it is 
readily confirmed that the turning points will move together 
as m increases (see Figure 6). It follows from (6-3) that if 


Ter +Ar, and r,-r,-Ar,, that Ag will decrease, an evident 


pattern in the table. 

It is noted that Table 3 was derived under the same 
conditions as the step index skew rays for Table 6 in Appendix 
E, except that here, a is not infinite but very big. The fact 
that the results of those two tables are approximately the 
same, demonstrates that the arbitrary a three dimensional 
analysis presented in this chapter is in agreement with the 


step index analysis. 
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D. FINAL REMARKS 

In this chapter, ray trajectories were established for the 
most general case of Chapter III. Although many numerical 
calculations have to be done, the model is still fairly 
simple. Also, the demonstration that trajectory dependence i 
the mode numbers agrees with theoretical expectations, is one 


more sign of the validity of the model. 
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VIt. QUASI ANALYTIC SOLUTION FOR SKEW RAYS WITH a=2 


A. INTRODUCTORY REMARKS 

In this chapter, the three-dimensional skew ray analysis 
will be examined for the special case a@=2. This case has 
special significance for telecommunications applications. It 
is for this profile, the so called parabolic profile, that 
multimode fibers exhibit an approximate minimum in modal 
dispersion [Ref. 1]. It turns out that the integrals 
describing the ray trajectories (see (6-1) and (6-2)) and 
radial mode number (see (3-15) with f(r)=0) can be evaluated 
analytically for the case a=2. Although numerical integration 
is no longer needed, the general relationship between m, «€(f), 
Y, and r, is too complex to be established without the 


assistance of a numerical root finding algorithm. 


B. ANALYSIS 


For a=2, (3-27) becomes 
(7-1) 


The radial mode number can be found by setting x-=r’ and 


f(r)=0 into (3-15), yielding 
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ee hes K-(2)2-¥ (7-2) 
m=— 7A) ea (=) x dx, 


where the turning points are (see (3-19)) given by 


<(B)| €7(B) -4( 7)? (753) 
X,2= - 








2 
After using [Ref. 18] 
[Vaxt+bxre Sf-- TEx sat +3) ef ax (7-4) 
as Vax?+bx+e xVax?+bx+c 
f dx 2 gine (7-5) 
Vax?+bx+c y-a b’-4ac 

and 
i aX ed gin (__bxr2e __ (7-6) 
xVax?+bx+c y-c |x| ¥b?- re 


in (7-2), it follows that 


€ (B) %-2(2)? : 
e@edn (es ee - 


-2X+e€ (B) 
a.) 7 
x| e?(B) -4 (5)? 


2 - V\2 
le (B) 4(—) 


sin! ( 





(7-7) 
where the limits are defined by (7-3). Equation (7-7) is not 
analytically invertible so it is necessary to have the 


computer find ¢(8) using the root searching algorithm shown in 
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Figure 19. Similarly, m,,, can be found after setting ¢(f8)=1 in 
(7-7) and directly evaluating the expression. 


In order to obtain the trajectory for a=2, (6-1) becomes 


[n?-NA? 
pets n,-NA € (B) E dx 
#(8) =" aa fy 


———————— (7-8) 
-& =. (252 
| X+E (B)X ia? 
which can be evaluated using (7-5) and leads to 
2 2 = 
— =, _ Yi -NA*e (8) Past -2X+e (B) E 
ia aoe Saami corer eres (7-9) 
2 = 2 
e*(B) (+) 
where xe[x,,X,]. The limits are given in (7-3). 
For $(x), (6-2) for a=2 becomes 
5d aL ga, ea 
ae) ae ie (7-10) 
%| -X+e (8B) X-(—)? 
Vv 
and by using (7-6), (7-10) becomes 
€(B)X-2(2)? 
o(®) =-5 [sin ES Sk (7-11) 


> 2 2 V\2 
x| 2) 4(+) 


where xe[x,,x,]. It is easy to show that substitution of (7-3) 


into (7-11) leads to 


$ (x,) =0 (7-12) 





o(%)=2, (7-13) 


respectively. After four applications of (7-11), $(x)=27, 
which proves the periodicity of ray trajectories with a=2, in 
agreement with References 8 and 20. 

Figure 28 shows a three-dimensional plot of a numerically 
calculated skew ray trajectory for conditions given on the 
Figure. Its axial view, Figure 29, demonstrates clearly that 
the ray trajectory on the XY plane is a closed orbit. As soon 
as the ray trajectory returns to the same point on the XY 
plane after four applications of (7-11), the z component of 
the ray trajectory is also forced to be periodic. 

Both the numerical ray tracing program used in Chapter VI, 
and the quasi analytic ray tracing program for a=2, were 
tested under the same conditions as given in Figure 28. It was 
determined that both programs produced exactly the same 
results with no greater then 0.1% discrepancy. These results 


lend confidence to the reliability of the models. 


C. EXAMPLES 

Table 5 was derived from repeated applications of (7-9) 
and (7-11) under optical fiber conditions consistent with V, 
NA, n, and a defined on Figure 28. The parameter Az is defined 
on Figure 28, and A¢ is defined on Figure 29. 

The discussion of the variation of parameters in Table 5 


has already been explained with similar tables in Chapter VI. 
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Figure 29: Axial view of Figure 28. 
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eee 
| 2 | 20 | so | 2 | 0.327 | 0.620 | 9.170 | 0.904 | 





Table 5: Variations of fiber parameters for a=2. 


It is easily checked that the parameter variations shown in 
the analytically generated Table 5 follow the same trends as 


seen in all the tables of Chapter VI. 


D. FINAL REMARKS 

In the previous chapter the case of a=2 was studied with 
the assistance of numerical integrations. It was demonstrated 
in Figure 27 that the trajectory actually followed a closed 
orbit for the specific fiber parameters. It was shown in this 
chapter that the result is exact and only depends on the 
condition that a=2. This somewhat remarkable result has been 
discussed in the literature. Specifically, the analogies 
between mechanics and optics provides the insight as the basis 


for this effect [Ref. 8]. It turns out that only two kinds of 
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potentials exist which can lead to closed orbits. First, there 


is the = potential, applicable to celestial bodies, and 


second, the r* potential associated with harmonic oscillators 


(Ref. 20]. It is the harmonic oscillator potential which is 


the analog of the parabolic profile. 











VIII. CONCLUSIONS 


The two main categories that deal with optical fibers are 
guided-wave electromagnetic modes and ray optics. Although 
these approaches both describe the propagation of waves, they 
are not directly connected. That is because the first category 
relies on waveguide theory to yield mode dispersion fields; 
while the second category uses initial conditions to establish 
ray trajectories. Previously, there was no way to move from 
angles of approach to waveguide modes, or relate ray 
trajectories to fields. The model presented in Chapter III has 
established a general formalism for connecting waveguide modes 
and ray trajectories. Specific cases are examined in Chapters 
IV through VII. 

As shown in all cases, the model relating the mode numbers 
is fairly simple. This and the fact that three-dimensional 
perspectives of ray trajectories can be easily produced, 
facilitates the generation of physical insight for a deeper 
understanding of optical fibers. 

It is noted that the results in Chapter V, Chapter VI and 
Appendix E show that the trajectory dependence on the mode 
numbers agrees with theoretical expectations. The direct proof 
in Appendix B that shows that the evolution equations are in 
agreement with standard Eikonal analysis makes the model 


presented valid. 
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Finally, the connection between modes and rays can 
ultimately lead into new methods in treating waveguide 
dispersion and bending losses, but this is a matter for 


further research. 
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APPENDIX A. CONCEPTUAL FLOWCHART 


The flowchart shown in Figure 30 summarizes the main steps 


in the development leading to mode specific trajectories. 





Input 
V, NA, a, n, 


| 
Predict 








Mmax>» ‘max 








Predict from m 


e(B), Fy Fo 








Ray trajectories are defined 
o(r), Z(r) 


Figure 30: Conceptual flowchart. 
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APPENDIX B. A COMPARISON TO EBIKONAL ANALYSIS 


The proposed evolution equations that describe the ray 
trajectories in the fiber optic according to the Eikonal 


approach (Ref. 9] are 





1 
ari mir) # 4)2 (B-1) 
az E r? 
and 
ap_ f (E-2) 
Zz r’ 
where E and ? are constants of the motion defined as 
E=-n(r) 2 (B-3) 
ds 
and 
ear? Sb (B-4) 


with s being the distance of the point measured along the ray 
trajectory. All other parameters appearing in (B-1) through 
(B-4), are defined in the discussion in Chapter III and Figure 
5. In the discussion to follow, we demonstrate the relation 
between the Eikonal constants of motion and the invariants £6 


and p. 


Definition of dl in Figure 5 and (3-4) leads to 





ds=dl. (B-5) 
By substituting (B-5) and (3-5) into (B-3), we have 
E=n(r)cos6, (B-6) 


which through combinations of (2-10) and (3-2) becomes 
B 
E==> B-7 
; ( ) 


which relates the phase constant $8 to the Eikonal constant E. 


Also, substitution of (3-5) and (3-7) into (B-4) gives 


g=y S8indsing (B-8) 
cos 


which when combined with (2-8), (2-10), (3-1) and (3-2) leads 


to 
( een (B-9) 


Equation (B-9) relates the azimuthal mode number pv and the 
phase constant 6 tc the Eikonal constant ?. 
Finally, substitution of (B-7) and (B-9) into (B-1) leads 
to 
Oe st 
dr : k,' (B-10) 
k?n?(r) -B?-(—)? 


and substitution of (B-9) and (B-10) into (B-2) in the form 


ar/dz 


ag/dz (B-11) 


leads to 
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adr or an SG (B-12) 
k’n? (zr) =6?=()? 
r 

Since (B-10) is equivalent to (3-28), and (B-12) is 
equivalent to (3-33), which are the basic evolution equations 
of the proposed analysis, there is an exact agreement between 
the analysis presented and the Eikonal analysis. Nevertheless, 
the treatment presented in this thesis connects Eikonal 


invariants to quantized wave mode numbers. 


75 





NOTATION 
Vv 
a 
ny 
ny 


NA 


kj 


k, 


€ (B) 


APPENDIX C. SYMBOL TABLE 
EXPLANATION 
V parameter 
core radius 
core index 
cladding index 
numerical aperture 


power profile parameter 


vacuum propagation constant (5) 


wn 
core propagation constant (—) 


w 
Cladding propagation constant ( 2) 
vacuum propagation wavelength (27/k) 

azimuthal mode number 


radial mode number 


coordinate system defining ray location (k,) 
angles which define ray orientation (k,) 


critical angle 
phase constant 


normalized phase constant 
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5 Sa SOME a EO extremes in radial excursions 
A trajectory period 
F(a) normalization integral 
2 2 
. ; ry - 
A index difference —! = 
2n; 
r,z ne 
aa 
“, E Eikonal invariants 
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APPENDIX D. PROPERTIES OF F(a) 


In this Appendix we examine F(a) given by (4-8), that is, 


F(a) ={,vi-u du. (D-1) 


As discussed in Chapter IV, F(a) cannot be solved analytically 
for all values of a. However, for some cases, solutions of 
F(a) can be obtained exactly. 


It is easy to show that 


F(0) =0 (D-2) 


and 
F(o) =1, (D-3) 


Also, by using [Ref. 18] 


[Vaxcrb ax= 21 (ox tee ‘ (D-4) 
a 


and 


x 2_ 4-2 2 
if a?-x? dx- NS + S sin! (2) ; (D-5) 


(D-1) takes the form 
F(1) =4 (D-6) 


and 
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_ 7 5 
F(Z) Ss (D-7) 


respectively. 


Since u%se"™, then for a<i 


2 3 
ut=1+aInu+ {einu)”, (olnu)” 


7 ..-=l+alnu, (D- 8) 


which when substituted into (D-1) leads to 


F, (a) ={ Yornu du. (D-9) 


After integration of (D-9) [Ref. 18], we find 


F (a) = Vom (D-10) 


for ol. 
For @ much greater than one, and since |u|<1, we may 


assume that 


y1-u* a1-3.. (D-11) 


Substitution of (D-11) into (D-1) and integration [Ref. 18] 


leads to 


2a+1 


D-12 
2a+2 


F,(a)= 





for a>l. 
Figure 8 shows a numerical calculation of F(a) for several 


values of a. Comparison with F,(q@) and F,(@) shows that the 
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functions are asymptotic to F(a) for small @ and large aq, 


respectively. 





APPENDIX E. 


REFERENCE 19 AS PRESENTED AT SSST 93 


A Heuristic Approach to the Computation of 3D-ray Trajectories 


in Step Index Optical Fibers 


Athanasios Nassopoulos and Ron Pieper 


Electrical and Computer Engineering 


Naval Postgraduate School 


Monterey, 


Abstract 


3D-ray trajectories in an optical fiber are 
derived through use of a simple correspondence 
principle. The ray trajectories are linked to the 
mode number of the exact waveguide solutions. 
The analysis and simulations presented will be in 
terms of dimensionless parameters needed to 
characterize the optical fiber. Specifically this 
the core index of 


includes the V-parameter, 


refraction, and the numerical aperture. The curves 
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CA 93943 


prepared are presented for allowed radial and 
azimuthal mode numbers and are presented in terms 
of dimensions which are normalized by the core 
radius. Although this approach is extendible to 
graded index fibers this development will not be 


presented in this report. 


introduction 


Exact waveguide solutions have been derived for z 


select number of index of refraction profiles. Among 





these exact solutions, the derivation for the 
step-index case commonly appears in most 
introductory and intermediate level optical fiber 
texts eg. [1]. Despite the mathematical 
importance of having an exact solution for the 
propagation in an optical fiber, the wave solutions 
can be difficult to interpret physically. It is not 
surprising that an alternate method to treating 
optical fibers, the approximate Eikonal ray 
approach [2], which in turn is based on Fermat's 
extremum principle, has been repeatedly applied in 
order to visualize effects which occur in optical 
fibers [3]. For example, bending losses [4} and 
dispersion effects [5] in optical fibers have been 
treated using this approach. However, in the 
Eikonal approach the derived ray trajectories are 
not automatically related to the radia! and 
azimuthal mode numbers. 

Recently an approach based on a simple 
correspondence between waves and rays has 
been proposed as a substitute for the first order 
WKB analysis in producing the Eikonal-equivalent 
of the ray propagation equations [6}. This 
correspondence principle has been used to 


qualitatively demonstrate that there is a direct 
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association between the exact wave modes and the 


allowed ray trajectories in an optical fiber. 


Background 


Starting from the Hemholtz equation which is the 
phasor form of the wave equation in linear, isotropic 


source-free homogeneous media, we have 


[v2+k2(r)]E-0, oa 


where k(f) is the wave number and E is the electric 


field. It is noted that this equation is only 


approximately valid in a graded index fiber. 
Consistent with the form given in (E-1) we assumed 


that 


k2()=n?(Hk? =k; +ky +z, (E-2) 

where k is the vacuum wave number, n(r} is the 
index of refraction, and r, @ and z are cylindrical 
coordinates needed to describe points in the fiber. 


Following a standard procedure [1] in waveguide 


analysis we assume that 


E{n=Rne'te hz (E-3) 


where v, the azimuthal mode number, is forced to 


be an integer due to periodic boundary conditions. 








of 


Figure 31: the 


problem 


Geometry 


B is the phase constant of the wave guide. Fir) is 
the radial solution to (E-1). In a step index fiber 
the exact waveguide solution for Fir) can be 
expressed on terms of Bessel Functions [1,3]. 
Using a simple correspondence rule, V-jk, 


it is easily shown [6] that 


(E-4) 


(E-5) 


A condition of phase synchronization [1] on the 
radial part of the wave vector results in an 


approximate integral expansion 


wt [2] K2p_p2-(¥)2 
rod pL 0-B (“ar 


where the radial mode number m takes on positive 


(E-6) 


integer values. 
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Equations (E-4), (E-5) and (E-6) become the basi: 
for the three dimensional ray modeling that will be 
the focus of the discussion that follows. It is notec 
that the components {k,, k,, k,) will describe thi 
direction of the ray and that (r, @, z) locate the ray 
For purposes of analysis the geometrica 
relationships are more carefully defined in the nex 


section. 


Geometry of the problem 


With reference to Fig. 31 the orientation of the 


wave vector is defined in terms of @ and &. | 


follows that 


k,=k(Asind sing (E-7) 
k,=K{)cosé (E-8) 
k,=k(Asin® cose. (E-9) 


Letting di represent the line variation in the ra 


location defined by 


: (E-10) 
di-dz 8,+dr'¥,+r db 8, 


and that the incremental components d@, dr, dz ai 


tangential to (E-7), (E-8) and (E-9) respectively, it follow 





dz=di cosé ({E-11) 
dr=al sin cost (E-12) 
f d=d/ sin® sing. (E-13) 


Through combinations of (E-7}-(E-13) the initial 


form of the trajectory equations becomes 


(E-14) 


kK tang (E-15) 
dr k, sak 
‘ k?n?(n -B°-(-? 


The evolution equations (E-14) and (E-15), that 
describe the ray trajectories, can be shown to be 
equivalent to the Eikonal solution [7]. Because of 
the 3D emphasis in this report we have chosen to 
concentrate on skew rays (v +0) instead of the 


less interesting meridional rays (v =O). 
Modeling Analysis. Step Index case 


We follow the usual step index notation 
nir)=n, for rsa and n(r)=n, for r>a.On Fig. 32 
the sum of the two terms in (E-6), k?n,?-(y/r)?, is 


shown plotted versus r. The condition that the 


(E-6) 


integration must remain positive would 
indicate that the hatched region in Fig. 32 defines 
the allowed extremities of the possible ray 


trajectories. It therefore follows that 


~1f* | p2n2 p22 
m fe | Ken p?-(~)Par. 


Specification of a radial number m and the 


(E-16) 


azimuthal mode numbers will dictate a specific 
phase constant §8,,, which are ordered [1] by 


convention 


ky22B yn? By mer>=-By mmax2 ke? (E-17) 


where the upper and lower bounds on £,,, are 
predictable from a combination of wave and ray 
analysis [1]. 

The turning point equation for r,,,, is obtained by 


setting the integrand in (E-16) to zero 


k?n,2-p2,,-(—-)#=0. (E-18) 
Fin 


The maximum allowed v is determined by setting 


min? et It 


Bum =k, and noting from Fig. 2 that r 


follows that 


Vnax=/nt [ka/n,?-n,7]=Int [(V| 


{E-19) 


where the quantity in the brackets is the well 
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Figure 32: Step index case 


known V parameter. The Int[ ] operator takes the 
integer part to produce a valid mode number.The 
general expression for r,, can be obtained from 


(E-18) as 


Fe ee : 
aah lash (E-20) 


where the second equation follows from geometric 
considerations. From (E-5) and (E-8), we find that 


(E-16) can be expressed as 


my)=lnt (Lf Ye neELy (E21) 
Tt Fran r 
where C=akn,sin@, r=r/a and  Crin=Cmin/a- 
Integration leads to 


my) =int (1 (/E7?-vsac (Sy). 2) 


85 


To find m,,,,.(¥) it is necessary to allow C to take its 
maximum value. As seen from (E-5) and (E-8) this 
occurs for the minimum 8. It follows from (E-17) 
that the maximum value of C is the V parameter. 


After direct substitution into (E-22), m,,.(v) 


satisfies 
Ming. (v) =Int (VE? vs00" ( My (E-23) 
where for vy =0 
Mpax(0) =i (E-24) 


1 
in agreement with a mode space analysis [8]. 
From the set of equations (E-20) and (E-22) it can 


easily be shown that 


mn =v | (1)2-1-sec( 1] (E-25) 
Troan Froin 


Given that m and p are known, the exact value of 


Tmin CaN be found by numerically solving (E-25)}. 


Normalized ray trajectories 


in this section the normalized trajectories for rays 
defined by a combination of standard optical 


parameters (V, NA, n,) and mode indices m and p. 








Figure 33: 3D representation showing 4 turning points 


Consistent with the purpose of this report the 
results for meridional rays (» =O) will be presented 
without significant discussion. 

Starting from (E-4) and (E-7) it easy to 


demonstrate that 


v=-f k,=-rk n, sin6 sing. (E-26) 


After substitution of (E-20) into (E-26) we have 


that 
(E-27) 


f, 
sing =--™, 
r 


which is needed in the analysis to follow. 
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The integral form of (E-15) which is 


sing dr 


@(A= f “tant dr= if : 


mv 1-shit, & 


(E-28) 


can be solved by substituting (E-27). Finally we 


have 


$()=sec(—L) {E-29) 
f, 


min 


which is not valid for »y =O. However from (E-15) 
and (E-27) g(r) =0 for »y =O. z(r) can be found by 


solving the integral form of(E-14) which is 





Figure 34: 
plane. 


Projection in polar 


ay=[ "22 or 


mn COSE 


(E-30) 


With a second application of(E-27) we find 


FG) =cney'F Pune 


where z=2/a. The angle @ can be easily calculated 


(E-31) 


from {E-20) as 


A, (E-32) 
‘min NA? 


and where r,,, can be obtained from (E-25). For 


min 


vy=0 it can be shown that 2(r)=fctn@ and 
sin@=(MNA)/(Vn,) are applicable to trajectory 
calculations of meridional rays. 

Along with the input parameters V, NA, n,, 
py (given v,,,, from (19) ) and m (given m,,,, from 
(E-23) ), the set of equations (E-25), (E-29), (E-31) 
and (E-32) the normalized 


determine ray 


trajectories of a skew ray in a step index fiber. 
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Discussion of results 


Figure 33 shows a 3D pilot of a numerically 
calculated skew ray trajectory for conditions given 
on the figure. The dense vertical lines are drawn at 


T=F,, and the low density lines are drawn at 


r=1.0. Points 1, 2, 3, 4 are turning points at 


T=1.0. Because of the 3D perspective the 
normalized increment in Az appears to change size. 
Figure 34 shows the corresponding axial view. The 
angle A@ is defined by rays joined at the turning 
points. 

Table 6 was derived from repeated applications 
under optical fiber conditions consistent with V, NA 
and n, defined on Figure 31. The trends relating 
variation of Ag and Az with »v and m are readily 
explained. tn the upper half of the table, for fixed m, 
the angle Ag@ increases with the _ rotational 
component, k,, as expected. Also, since m (E-6) is 
constant, # must decrease. This is manifest by AZ 
decreasing with v. In the lower half of the table, 
where v is fixed, as m increases, 8, as seen from 


(E-6) must again drop. Again AZ must therefore 


decrease. And finally, the parallel decrease in Ag (E- 


17.25 
12.57 
09.84 
12.37 
07.38 
06.10 


Table 6:Representative variation 


in ray parameters with wave mode 


numbers. 

29) can be explained geometrically as seen from 
figure 34 by a required decrease in r,,,. This is 
apparent by noting that as the radial mode number 
increases with pv fixed the phase constant 8 
should, as seen from (E-6) decrease. Conditions 
(E-5) and (E-8) then indicate @ will increase and 


therefore from (E-20) r,,. will decrease. 


min 


Conclusions 


The main point of this report on optical 
fiber ray calculations is to demonstrate the 
existence of a pedagogically attractive alternative 


to more formal methods such as WKB or Eikonal 


analysis. An important feature of the method 
presented is that it links the traditional cylindrical 
waveguide modenumbers for the exact solution to 
the ray trajectories. Because of this link the 
numerical implementation of the model can provide 
valuable physical insight for the fiber optic 
application areas such as sensors or lightwave 


telecommunications. 
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